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Abstract 

We describe in detail the theory underpinning the measurement of density matrices of a pair of quantum 
two- level systems ("qubits"). Our particular emphasis is on qubits realized by the two polarization de- 
grees of freedom of a pair of entangled photons generated in a down-conversion experiment; however the 
discussion applies in general, regardless of the actual physical realization. Two techniques are discussed, 
namely a tomographic reconstruction (in which the density matrix is linearly related to a set of measured 
quantities) and a maximum likelihood technique which requires numerical optimization (but has the ad- 
vantage of producing density matrices which are always non-negative definite) . In addition a detailed error 
analysis is presented, allowing errors in quantities derived from the density matrix, such as the entropy or 
entanglement of formation, to be estimated. Examples based on down-conversion experiments are used to 
illustrate our results. 
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I. INTRODUCTION 



The ability to create, manipulate and characterize quantum states is becoming an increasingly important area 
of physical research, with implications for areas of technology such as quantum computing, quantum cryptography 
and communications. With a series of measurements on a large enough number of identically prepared copies of an 
quantum system, one can infer, to a reasonable approximation, the quantum state of the system. Arguably, the first 
such experimental technique for determining the state of quantum system was devised by George Stokes in 1852 [Q. 
His famous four parameters allow an experimenter to determine uniquely the polarization state of a light beam. With 
the insight provided by nearly 150 years of progress in optical physics, we can consider coherent light beams to be an 
ensemble of two-level quantum mechanical systems, the two levels being the two polarization degrees of freedom of the 
photons; the Stokes parameters allow one to determine the density matrix describing this ensemble. More recently, 
experimental techniques for the measurement of the more subtle quantum properties of light have been the subject of 
intensive investigation (see ref. for a comprehensive and erudite exposition of this subject). In various experimental 
circumstances is has been found reasonably straightforward to devise a simple linear tomographic technique in which 
the density matrix (or Wigner function) of a quantum state is found from a linear transformation of experimental 
data. However, there is one important drawback to this method, in that the recovered state might not correspond 
to a physical state because of experimental noise. For example, density matrices for any quantum state must be 
Hermitian, positive semi-definite matrices with unit trace. The tomographically measured matrices often fail to be 
positive semi-definite, especially when measuring low-entropy states. To avoid this problem the "maximum likelihood" 
tomographic approach to the estimation of quantum states has been developed (|-|7). In this approach the density 
matrix that is "mostly likely" to have produced a measured data set is determined by numerical optimization. 

In the past decade several groups have successfully employed tomographic techniques for the measurement of 
quantum mechanical systems. In 1990 Risley et al. at North Carolina State University reported the measurement 
of the density matrix for the nine sublevels of the n = 3 level of hydrogen atoms formed following collision between 
H"*" ions and He atoms, in conditions of high symmetry which simplified the tomographic problem [^. Since then, 
in 1993 Smithey et al. (University of Oregon), made a homodyne measurement of the Wigner function of a single 
mode of light |^ . Other explorations of the quantum states of single mode light fields have been made by Mlynek et 
al. (University of Konstanz, Germany) and Bachor et al. (Australian National University) jl^. Other quantum 
systems whose density matrices have been investigated experimentally include the vibrations of molecules |]l2|, the 
motion ions and atoms p^ , p^ and the internal angular momentum quantum state of the F = 4 ground state of a 
cesium atom [ p^ . The quantum states of multiple spin- half nuclei have been measured in the high-temperature regime 
using NMR techniques jl^] , albeit in systems of such high entropy that the creation of entangled states is necessarily 
precluded p^ . The measurement of the quantum state of entangled qubit pairs, realized using the polarization degrees 
of freedom of a pair of photons created in a parametric down-conversion experiment was reported by us recently [ p^ . 

In this paper we will examine techniques in detail for quantum state measurement as it applies to multiple correlated 
two- level quantum mechanical systems (or "qubits" in the terminology of quantum information). Our particular 
emphasis is qubits realized via the two polarization degrees of freedom of photons, data from which we use to 
illustrate our results. However, these techniques are readily applicable to other technologies proposed for creating 
entangled states of pairs of two-level systems. Because of the central importance of qubit systems to the emergent 
discipline of quantum computation, a thorough explanation of the techniques needed to characterize the qubit states 
will be of relevance to workers in the various diverse experimental fields currently under consideration for quantum 
computation technology ||l^. This paper is organized as follows: In Section II we explore the analogy with the Stokes 
parameters, and how they lead naturally to a scheme for measurement of an arbitrary number of two-level systems. 
In Section HI, we discuss the measurement of a pair of qubits in more detail, presenting the validity condition for an 
arbitrary measurement scheme and introducing the set of 16 measurements employed in our experiments. Section IV 
deals with our method for maximum likelihood reconstruction and in Section V we demonstrate how to calculate the 
errors in such measurements, and how these errors propagate to quantities calculated from the density matrix. 



II. THE STOKES PARAMETERS AND QUANTUM STATE TOMOGRAPHY 

As mentioned above, there is a direct analogy between the measurement of the polarization state of a light beam 
and the measurement of the density matrix of an ensemble of two-level quantum mechanical systems. Here we explore 
this analogy in more detail. 
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A. Single qubit tomography 



The Stokes parameters are defined from a set of four intensity measurements [g^ : (i) witli a filter that transmits 
50% of the incident radiation, regardless of its polarization; (ii) with a polarizer that transmits only horizontally 
polarized light; (iii) with a polarizer that transmits only light polarized at 45° to the horizontal; and (iv) with a 
polarizer that transmits only right circularly polarized light. The number of photons counted by a detector, which is 
proportional to the classical intensity, in these four experiments are as follows: 

no = y mm) + {v\p\v)) = y mm) + {Lm)) 
n,=Mmm)) 
= Y mm) + {Rm) + mm) + mm)) 

n2=N{{Dm)) 

= y mm) + mm) - imm) + mm)) 
mm)) 



(2.1) 



Here \H), \V), \D) = (|i7) - \V)) /^/2 = exp(i7r/4) (|i?) + i\L)) /V2 and \R) = (|i/) - i\V)) /V2 are the kets rep- 
resenting qubits polarized in the linear horizontal, linear vertical, linear diagonal (45°) and right-circular senses 
respectively, p is the (2 x 2) density matrix for the polarization degrees of the light (or, for a two-level quantum 
system) and Af is constant dependent on the detector efficiency and light intensity. The Stokes parameters, which 
fully characterize the polarization state of the light, are then defined by 
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(2.2) 



We can now relate the Stokes parameters to the density matrix p by the formula 



1 Si ~ 
^7^0^" 



(2.3) 



where cto = \R){R\ + \L){L\ is the single qubit identity operator and cti = \B){L\ + \L){R\, a2 = i{\L) {R\ — \R) {L\ and 
tT3 = \R){R\ — \L){L\ are the Pauli spin operators. Thus the measurement of the Stokes parameters can be considered 
equivalent to a tomographic measurement of the density matrix of an ensemble of single qubits. 



B. Multiple beam Stokes parameters: multiple qubit tomography 

The generalization of the Stokes scheme to measure the state of multiple photon beams (or multiple qubits) is 
reasonably straightforward. One should, however, be aware that importance differences exist between the one-photon 
and the multiple photon cases. Single photons, at least in the current context, can be described in a purely classical 
manner, and the density matrix can be related to the purely classical concept of the coherency matrix [^T| . For mul- 
tiple photon one has the possibility of non-classical correlations occurring, with quintessentially quantum-mechanical 
phenomena such as entanglement being present. We will return to the concept of entanglement and how it may be 
measured later in this paper. 

An n-qubit state is characterized by a density matrix which may be written as follows: 

1 ^ 

^^=2" ^ ''»i,*2,...»„o-,^ (g)(T,2 (g) . . . «) o-,„, (2.4) 

ii,i2,---in—0 
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where the 4" parameters ri-^^i^^,,,i^ are real numbers. The normaUzation property of the density matrices requires that 
''0,0,... = 1, and so the density matrix is specified by 4" — 1 real parameters. The symbol ® represents the tensor 
product between operators acting on the Hilbert spaces associated with the separate qubits. 

As Stokes showed, the state of a single qubit can be determined by taking a set of four projection measurements 
which are represented by the four operators /io = \H){H\ + \V){V\, p-i = \H){H\, (12 = \D){D\, fi^ = \R){R\. 
Similarly, the state of two qubits can be determined by the set of 16 measurements represented by the operators 
Ai ® Aj ihj = 0,1,2,3). More generally the state of an n-qubit system can be determined by 4" measurements 
given by the operators Aii ® Ai2 ^ • • ■ ^ Ai„ = 0; li 2, 3) and (fc = 1, 2, . . . n). This 'tree' structure for multi-qubit 
measurement is illustrated in fig.l. 

The proof of this conjecture is reasonably straightforward. The outcome of a measurement is given by the formula 



= AATr{pA}: 



(2.5) 



where p is the density matrix, jl is the measurement operator and A/" is a constant of proportionality which can be 
determined from the data. Thus in our n-qubit case, the outcomes of the various measurement are 



A/'Tr {p {jl,, ® Pi. 



(2.6) 



Substituting from eq.(2.4) we obtain 



2" 



ii,i2,...jn=o 



(2.7) 



As can be easily verified, the single qubit measurement operators Ai are linear combinations of the Pauli operators 
aj, i.e. Pi = X^'i^o ^ij^j' where Tij are the elements of the matrix 



/I \ 

1/2 1/2 
1/2 1/2 
\ 1/2 1/2 J 

Further, we have the relation Ti {ai<7j} — 26ij (where 6ij is the Kronecker delta). Hence eq.(^) becomes 
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.T 



Introducing the left-inverse of the matrix T, defined so that X^fc^oC^ ^)ik^kj — 5ij and whose elements are 



/ 1 0\ 

-12 
-10 2 
\ -1 2 / 



we can find a formula for the parameters ri^^i^^,,,i^ in terms of the measured quantities 



''Zl,22,...2/i 7 



(2.8) 



(2.9) 



(2.10) 



= .i'y -i^ • 



(2.11) 



In eq.( ^.ll ) we have introduced the n-photon Stokes parameter Sij^^i2,...i„, defined in an analogous manner to the 
single photon Stokes parameters give in eq.(2.2). 

Since, as already noted, ro,o....o = 1, we can make the identification iSo,o,...o = A/", and so the density matrix for the 
n-qubit system can be written in terms of the Stokes parameters as follows: 



On / ^ 



5,- 



2" . . ^ 5o,o....o 



(2.12) 
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This is a recipe for measurement of the density matrices which, assuming perfect experimental conditions and the 
complete absence of noise, will always work. It is important to realize that the set of four Stokes measurements 
{/to, Aij A2j As} &re not unique: there may be circumstances in which it is more convenient to use some other set, 
which are equivalent. A more typical set, at least in optical experiments, is Ao = \H){H\, fi[ = |V^)(y|, /ij = \D){D\, 
fi',^\R){R\. 

In the following section we will explore more general schemes for the measurement of two qubits, starting with a 
discussion, in some detail, of how the measurements are actually performed. 



III. GENERALIZED TOMOGRAPHIC RECONSTRUCTION OF THE POLARIZATION STATE OF TWO 

PHOTONS 



A. Experimental set-up 

The experimental arrangement used in our experiments is shown schematically in Fig.l. An optical system consisting 
of lasers, polarization elements and non-linear optical crystals (and collectively characterized for the purposes of this 
paper as a "black-box",) is used to generate pairs of qubits in an almost arbitrary quantum state of their polarization 
degrees of freedom. A full description of this optical system and how such quantum states can be prepared can be 
found in ref. [p2|-p4| |^ The output of the black box consists of a pair of beams of light, whose quanta can be measured 
by means of photo-detectors. To project the light beams onto a polarization state of the experimenter's choosing, 
three optical elements are placed in the beam in front of each detector: a polarizer (which transmits only vertically 
polarized light), a quarter- wave plate and a half- wave plate. The angles of the fast axes of both of the waveplates can 
be set arbitrarily, allowing the \V) projection state fixed by the polarizer to be rotated into any polarization state 
that the experimenter may wish. 

Using the Jones calculus notation, with the following convention, 

J)=|V), (J)=|H), (3.1) 

where |V) (|H)) is the ket for a vertically (horizontally) polarized beam, the effect of quarter and half wave plates 
whose fast axes are at angles q and h with respect to the vertical axis, respectively, are given by the 2x2 matrices 

fj 1 f i- C0s(2q) sin(2g) 

UQWP{q}-^^ sin(2<z) z + cos(2g) 

fr _ f cos(2/i) - sin(2/i) \ , , 

UHWP{q) - _ ^i^(2/i) - cos(2/.) ) ■ ^^-2) 

Thus the projection state for the measurement in one of the beams is given by 

\4'roAh,q)) - UQWp{q) ■ UHWp{h) ■ \ 

^a{h,qm+b{h,q)\V), (3.3) 
where, neglecting an overall phase, the functions a(ft,, q) and 6(/i, q) are given by 

a{h, q) — — ^ (sin(2/i) — i sin[2(/i — q)]) , 
v2 

b{h, g) = — ^ (cos(2/i) + i cos[2(/i - q)]) . (3.4) 
The projection state for the two beams is given by 



^It is important to realize that the entangled photon pairs are produced in a non- deterministic manner: one cannot specify 
with certainly when a photon pair will be emitted; indeed there is a small probability of generating four, or six or higher 
numbers of photons. Thus we can only post-selectively generate entangled photon pairs: i.e. one only knows that the state was 
created after if has been measured. 
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= a(/ii,gi)a(/i2,g2)|HH) +a(/ti,gi)6(/i2,92)|HV) + 

6(/ii,gi)a(/i2,g2)|VH) +&(/ii,gi)&(/i2,g2)|VV). (3.5) 

We shall denote the projection state corresponding to one particular set of waveplate angles {/ii,!^, 91,2^, ft.2.1/, '?2.i/} 
^ by the ket thus the projection measurement is represented by the operator fi^ = \'ip^){'ipi^\. Consequently, the 
average number of coincidence counts that will be observed in a given experimental run is 

=AA(V',|p|^/.,) (3.6) 

where p is the density matrix describing the ensemble of qubits, and A/" is a constant dependent on the photon flux 
and detector efficiencies. In what follows, it will be convenient to consider the quantities Si, defined by 

s„ ^ {ipi.\p\ip„)- (3.7) 



B. Tomographically Complete set of Measurements 

In Section II we have given one possible set of projection measurements which which uniquely determine 

the density matrix p. However, one can conceive of situations in which these will not be the most convenient set of 
measurements to make. Here we address the problem of finding other sets of suitable measurements. The smallest 
number of states required for such measurements can be found by a simple argument: there are 15 real unknown 
parameters which determine a 4 x 4 density matrix, plus there is the single unknown real parameter A/", making a 
total of 16. 

In order to proceed it is helpful to convert the 4x4 matrix p into a 16-dimensional column vector. To do this we 
use a set of 16 linearly independent 4x4 matrices {T^} which have the following mathematical properties: 



Tr ■ ^m} = ^i^.fj^ 

16 

i = ^f,Tr{f, -i} Vi, (3.8) 

1/=! 

where A is an arbitrary 4x4 matrix. Finding a set of T^, matrices is in fact reasonably straightforward: for example, 
the set of (appropriately normalized) generators of the Lie algebra SU{2) (g) SU{2) fulfill the required criteria (for 
reference, we list this set in Appendix A). These matrices are of course simply a re-labeling of the two-qubit Pauli 
matrices ai (g) d'j ~ 0, 1, 2, 3) discussed above. Using these matrices the density matrix can be written as 

16 

p = Y.^.r,, (3.9) 

1^=1 

where ri, is the v-th element of a sixteen element column vector, given by the formula 



Tr 



{f . • p} (3.10) 



Substituting from eq.( |3.9D into eq.(3.6), we obtain the following linear relationship between the measured coincidence 
counts and the elements of the vector r^: 

16 

=AA^B,,^r^ (3.11) 



"•"Here the first subscript on the waveplate angle refers one of the two photon beams; the second subscript distinguishes which 
of the sixteen different experimental states is under consideration. 
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where the 16 x 16 matrix -Bjy,^ is given by 



(3.12) 



Immediately we find a necessary and sufiicient condition for the completeness of the set of tomographic states {IV"!^)}: 
if the matrix is non- singular, then eg. ( 3. 11 ) can be inverted to give 



(3.13) 



The set of sixteen tomographic states which we employed are given in Table 1. They can be shown to satisfy the 
condition that B^^^ is non- singular. By no means are these states unique in this regard: these were the states chosen 
principally for experimental convenience. 

These states can be realized by setting specific values of the half- and quarter- wave plate angles. The appropriate 
values of these angles (measured from the vertical) are given in Table 1. Note that overall phase factors do not affect 
the results of projection measurements. 

Substituting eq.(3.13) into eq.(3.9), we find that 

16 

/5= (AA)-i^M,n, 
1^=1 

16 

= Y^aXs,, (3.14) 



where the sixteen 4x4 matrices Mi, are defined by 



16 



M, 



(3.15) 



v=l 



The introduction of the M^, matrices allows a compact form of linear tomographic reconstruction, eq.( |3.14 ), will be 
most useful in the error analysis that follows. These Mi, matrices, valid for our set of tomographic states, are listed 
in Appendix B, together with some of their important properties. We can use one of these properties, eq.(B.6), to 
obtain the value of the unknown quantity Af. That relationship implies 



^Tr{M,}|^,)(^,|p = p. 
Taking the trace of this formula, and multiplying by M we obtain: 

For our set of tomographic states, it can be shown that 

^Tr{Af.}: 

hence the value of the unknown parameter M in our experiments is given by: 



1 if J/ =1,2, 3, 4 
if i/ = 5,...16, 



(3.16) 



(3.17) 



(3.18) 



TV = 

= N{{HH\p\HH) + {HV\p\HV) + {VH\p\VH) + {VV\p\VV)) . 
Thus we obtain the final formula for the tomographic reconstruction of the density matrices of our states: 



(3.19) 



(3.20) 
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As an example, the following set of 16 counts were taken for the purpose of tomographically determining the density 
matrix for an ensemble of qubits all prepared in a specific quantum state: ni — 34749, 712 = 324,^3 — 35805,714 — 
U4,n5 = 16324,716 = 17521, 717 = 13441,778 = 16901, Tig = 17932, 77io = 32028, 77ii = 15132, 7712 = 17238, 7113 === 



13171,7714 = 17170,7715 = 16722,7716 = 33586 Applying eq.(3.20) we find 



/ 0.4872 -0.0042 + ^0.0114 -0.0098 - i0.0178 0.5192 + ^0.0380 \ 

-0.0042 -iO.0114 0.0045 0.0271-^0.0146 -0.0648 - i0.0076 

-0.0098 + ^0.0178 0.0271 + ^0.0146 0.0062 -0.0695 + 70.0134 

\ 0.5192-^0.0380 -0.0648 + ^0.0076 -0.0695 - i0.0134 0.5020 / 



(3.21) 



This matrix is shown graphically in fig. 3a. 

Note that, by construction, the density matrix is normalized, i.e. Tr{/5} — 1 and Hermitian, i.e. — p . However, 
when one calculates the eigenvalues of this measured density matrix, one finds the values 1.02155, 0.0681238, - 
0.065274 and -0.024396; and also, Tr{p^} = 1.053 . Density matrices for all physical states must have the property 
of positive semi-definiteness, which (in conjunction with the normalization and Hermiticity properties) imply that 
all of the eigenvalues must lie in the interval [0,1], their sum being 1; this in turn implies that < Tr{/5^} < 1. 
Clearly the density matrix reconstructed above by linear tomography violates these condition. From our experience 
of tomographic measurements of various mixed and entangled states prepared experimentally, this seems to happen 
roughly 75% of the time for low entropy, highly entangled states; it seems to have a higher probability of producing 
the correct result for states of higher entropy, but the cautious experimenter should check every time. The obvious 
culprit for this problem is experimental inaccuracies and statistical fluctuations of coincidence counts, which mea n 
that the actual numbers of counts recorded in a real experiment differ from those that can be calculated by eq. ( t3.6[ ) . 
Thus the linear reconstruction is of limited value for states of low entropy (which are of most experimental interest 
because of their application to quantum information technology); however as we shall see, the linear approach does 
provide a useful starting point for the numerical optimization approach to density matrix estimation which we will 
discuss in the next section. 

IV. MAXIMUM LIKELIHOOD ESTIMATION 

As mentioned in Section III, the tomographic measurement of density matrices can produce results which violate 
important basic properties such as positivity. To avoid this problem, the maximum likelihood estimation of density 
matrices may be employed. Here we describe a simple realization of this technique. 

A. Basic approach 

Our approach to the maximum likelihood estimation of the density matrix is as follows: 

(i) Generate a formula for an explicitly "physical" density matrix, i.e. a matrix which has the three important 
properties of normalization, Hermiticity and positivity. This matrix will be a function of 16 real variables (denoted 
{ti,t2, • ■ -^le})- We will denote the matrix as Pp(ti,t2, ■ ■ - tie). 

(ii) Introduce a "likelihood function" which quantifies how good the density matrix Pp{ti, ^2, • • • ^le) is in relation to 
the experimental data. This likelihood function is a function of the 16 real parameters t^, and of the 16 experimental 
data n^. We will denote this function as C{ti,t2, ■ . ■ tie; ni,n2, ■ ■ ■ uiq). 

(iii) Using standard numerical optimization techniques, find the optimum set of variables {t'f^*\t'^^*\ . . . t^iQ*^} for 
which the function C{ti, ^2, • ■ ■ ^le! "-ij n-2, ■ ■ ■ nie) has its maximum value. The best estimate for the density matrix is 
then p(4°P*\4°P*\...4f )). 

The details of how these three steps can be carried out are described in the next three sub-sections. 

B. Physical Density Matrices 

The property of non-negative definiteness for any matrix Q is written mathematically as 

{ip\g\i^)>o viv^). (4.1) 
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Any matrix that can be written in the form Q — T^T must be non-negative definite. To see that this is the case, 
substitute into eq.([4.l|): 



(4.2) 



where we have defined = T'lV')- Furthermore (f tr)t = f t(f t)t = f ^T, i.e. = f tf must be Hermitian. To 
ensure normahzation, one can simply divide by the trace: thus the matrix g given by the formula 
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= T^f /Tr{f tf } 



(4.3) 



has all three of the mathematical properties which we require for density matrices. 

For the two qubit system, we have a 4 x 4 density matrix with 15 independent real parameters. Since it will be 
useful to be able to invert relation ( |4.3| ), it is convenient to choose a tri-diagonal form for T: 



( 


h 







h 


+ ite h 








hi 


+ iti2 h + its 









+ iti6 ti3 + itii 


ig + itiQ 


ti j 



fit) 



Thus the explicitly "physical" density matrix pp is given by the formula 

Ppit)^fHt)m/TT{fHt)fit)}. 



(4.4) 



(4.5) 



For future reference, the inverse relationship, by which the elements of T can be expressed in terms of the elements 
of p, is as follows: 



/ 



T ^ 











\ 



Mil], 



M 



(2) 
11,23 



/P44 



P41 

VP44 



/p^V^il 

P42 
VP44 



P44 



(4.6) 



/P44 



fpii J 



Here we have used the notation A — Det(/5); Mfj' is the first minor of p, i.e. the determinant of the 3x3 matrix 

formed by deleting the i-th row and j-th column of p; M\j is the second minor of p, i.e. the determinant of the 
2x2 matrix formed by deleting the i-th and fc-th rows and j-th and l-ih columns oi p {i ^ k and j ^ I)- 



C. The Likelihood Function 



The measurement data consists of a set of 16 coincidence counts Ui, {v = 1,2,... 16) whose expected value is 
fiy ~ M{i]j^\p\%pi,). Let us assume that the noise on these coincidence measurements has a Gaussian probability 
distribution. Thus the probability of obtaining a set of 16 counts {rii, ^2, . . . A^ie} is 



P{ni,n2, ■ ■ - nie) = 



1 



16 



Norm 



n "^^p 



2ol 



(4.7) 



where is the standard deviation for i^-th coincidence measurement (given approximately by \fn^) and Norm is the 
normalization constant. For our physical density matrix pp the number of counts expected for the v-th measurement 
is 

(ii, i2, ■ ■ • tie) = N{iJu\Pp {ti,t2, . . . tie) li'^)- (4.8) 
Thus the likelihood that the matrix pp {ti,t2, ■ ■ ■ tie) could produce the measured data {rii, n2, . . . nie} is 
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P(ti,t2, ■ --tie) = ^^^^ Yl "^^P 

1^=1 



'' {M{')pi,\pp (ti, t2, . . . ^le) lipu) - n^,) 



2Af{ipi,\pp {ti,t2, ■ ..he) lip^) - n^) 



(4.9) 



where M = X^^^i ^v- 

Rather than find maximum value of P (ti, t2, . . . tig) it simplifies things somewhat to find the maximum of its 
logarithm (which is mathematically equivalent) ^ Thus the optimization problem reduces to finding the minimum of 
the following function: 

r(. . . ^ {■N'{^v\Pp{tiM. ■ ■ ■tw)\'>l^u) -n^f ,„^ 

This is the "likelihood" function which we employed in our numerical optimization routine. 

D. Numerical Optimization 

We used the Mathematica 4.0 routine FindMinimum which executes a multidimensional Powell direction set algo- 
rithm (see ref. | p5[ | for a description of this algorithm). To execute this routine, one requires an initial estimate for 
the values of , i2 , ■ ■ • ^le • For this, we used the tomographic estimate of the density matrix in the inverse relation 



(4.6), allowing us to determine a set of values for ^1,^2, • • -iie- Since the tomographic density matrix may not be 
non-negative definite, the values of the t^'s deduced in this manner are not necessarily real. Thus for our initial guess 
we used the real parts of the t^^s deduced from the tomographic density matrix. 
For the example given in Section 2, the maximum likelihood estimate is 

1 0.4833 -I- i0.0329 \ 
_n noQR _ ,n nn77 

(4.11) 



( 0.5069 -0.0239 -I- i0.0106 -0.0412 - i0.0221 0.4833 -I- i0.0329 \ 

-0.0239 -i0.0106 0.0048 0.0023 -I- i0.0019 -0.0296 - i0.0077 

-0.0412 -I- i0.0221 0.0023 -i0.0019 0.0045 -0.0425 -I- i0.0192 

V 0.4833 - i0.0329 -0.0296 -I- i0.0077 -0.0425 - i0.0192 0.4839 / 



This matrix is illustrated in Fig. 2b. In this case, the matrix has eigenvalues 0.986022, 0.0139777, and 0; and 
Tr{/5^} — 0.972435, indicating, while the linear reconstruction gave a non-physical density matrix, the maximum- 
likelihood reconstruction gives a legitimate density matrix. 

V. ERROR ANALYSIS 

In this section we present an analysis of the errors inherent in the tomographic scheme described in Section III. 
Two sources of errors are found to be important: the shot noise error in the measured coincidence counts and the 
uncertainty in the settings of the angles of the waveplates used to make the tomographic projection states. We will 
analyze these two sources separately. 

In addition to determining the density matrix of a pair of qubits, one is often also interested in quantities derived 
from the density matrix, such as the entropy or the entanglement of formation. For completeness, we will also derive 
the errors in some of these quantities. 

A. Errors due to Count Statistics 

From eq. (|3.20| ) we see that the density matrix is specified by a set of sixteen parameters Si, defined by 

Sy^n^/N, (5.1) 



^Note that here we neglect the dependence of the normalization constant on t\,t2, ■ ■ ■ tie, which only weakly effects solution 
for the most likely state 
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where rij, are the measured coincidence counts and JV — X]i/=i'^i'- We can determine the errors in using the 
foUowing formula Efl] 



16 

E 

A.K=1 



driK, 



SnxSn,^, 



(5.2) 



where the over-bar denotes the ensemble average of the random uncertainties Ss^, and Sn\. The measured coincidence 
counts nx are statistically independent Poissonian random variables, which implies the following relation: 



SnxSn^ = nxSx,K, 



where 6x « is the Kronecker delta. 



Taking the derivative of eq.(5.1), we find that 

ds 



where 



A=l 



1 if 1 < 1/ < 4 
if 5 < < 16. 



(5.3) 



(5.4) 



(5.5) 



Substituting from eq.(5.4) into eq.(5.2) and using eq.(5.3), we obtain the result 



(5.6) 



In most experimental circumstances Af ^ 1, and so the second term in eq.( |5.6| ) is negligibly small in comparison to 
the first. We shall therefore ignore it, and use the approximate expression in the subsequent discussion; 



OSvOS^ ~ -JJ^Ou^fi = -^Oiy.p. 



(5.7) 



B. Errors due to Angular Settings Uncertainties 

Using the formula (^^) for the parameters Su we can find the dependence of the measured density matrix on errors 
in the tomographic states. The derivative of Su with respect to some generic wave-plate settings angle 9 is 



de 



where is the kct of the u-th projection state [see eq.(3^)]. Substituting from eq.( ^3.14| ) wc find 



89 



89 



(5.8) 



(5.9) 



For convenience, we shall label the four waveplate angles {hi ,y,qi^,h2.^,q2,u} which specify the v-th state by 
{^'i/.i, ^'i/,2, ^'i/,3, 6*1^,4} respectively. Clearly the /i-th state does not depend on any of the v-th set of angles. Thus we 
obtain the following expression for the derivatives of with respect to waveplate settings: 



8s 



16 



89 



— "J'.A ^fiJ L',fJ.^ 

M=l 



(5.10) 



where 



89,, 



89,. 



(5.11) 
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The 1024 quantities f^^l can be determined by taking the derivatives of the functional forms of the tomographic states 



given by eqs.(3^) and (^^), and evaluating those derivatives at the appropriate values of the arguments (see Table 
l). ^ 

The errors in the angles are assumed to be uncorrelated, as would be the case if each wave-plate were adjusted for 
each of the 16 measurements. In reality, for qubit experiments, only one or two of the four waveplates are adjusted 
between one measurement and the next. However the assumption of uncorrelated angular errors greatly simplifies the 
calculation (which is, after all, only an estimate of the errors), and seems to produce reasonable figures for our error 
bars Q Thus with the following assumption 

M~M~^s,,^s,^,iAef, (5.12) 

(where A9 is the RMS uncertainty in the setting of the waveplate, with an estimated value of 0.25° for our apparatus) 
we obtain the following expression for the errors in due to angular settings: 



16 



= 1 e,A=l 



Combining eqs. (|l3|) and we obtain the following formula for the total error in the quantities s^: 



where 



A. 



i=\ £,A=1 



(5.13) 



(5.14) 



(5.15) 



These sixteen quantities can be calculated using the parameters Si, and the constants /i^i. Note that the same result 
can be obtained by assuming a 'priori that the errors in the are all uncorrelated, with = 8s\\ the more rigorous 
treatment given here is however necessary to demonstrate this fact. For a typical number of counts, say M = ooffff it 
is found that the contribution of errors from the two causes is roughly comparable; for larger numbers of counts, the 
angular settings will become the dominant source of error. 

Based on these results, the errors in the values of the various elements of the density matrix estimated by the linear 
tomographic technique described in Section 3 are as follows: 



^ ds^ as,; 



iy,^— 1 
16 

E(*^. 



■A. 



(5.16) 



where M^(i j-^ is the i, j element of the matrix M^,. 

A convenient way in which to estimate errors for a maximum likelihood tomographic technique (rather than a linear 
tomographic technique) is to employ the above formulae, with the slight modification that the parameter should 
be recalculated from eq.(3.7) using the estimated density matrix pest- This does not take into account errors inherent 
in the maximum likelihood technique itself. 



C. Errors in Quantities Derived from the Density Matrix 



When calculating the propagation of errors, it is actually more convenient to use the errors in the s^, parame- 
ters, (given by eq.(5.15), rather than the errors in the elements of density matrix itself (which have non-negligible 
correlations). 



**In other experimental circumstances, such as the measurement of the joint state of two spin 1/2 particles, the tomography 
would be realized by performing unitary operations on the spins prior to measurement. In this case, an assumption analogous 
to ours will be wholly justified. 
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1. von Neumann Entropy 

The von Neumann entropy is an important measure of the purity of a quantum state p. It is defined by 

5 = -Tr{plog2 m 

4 

^-^Pa^Og^Pa, (5.17) 

a=l 

where Pa is an eigenvalue of p, i.e. 

p\q^a) ^Pa\(ba), (5.18) 

1 0a) being the a-th eigenstate (a= 1,. . .4). The error in this quantity is given by 

16 .^^.2 



v=l 



dS \ / dpa \ f dS 



Applying the chain rule, we find 

;;:)-i:(i:j(£j. 

The partial differential of a eigenvalue can be easily found by perturbation theory. As is well known (e.g. |Q) the 
change in the eigenvalue Xa of a matrix W due to a perturbation in the matrix SW is 

6Xa = {(j)a\6W\q}a), (5.21) 

where {(fia) is the eigenvector of W corresponding to the eigenvalue Xa- Thus the derivative of Xa with respect to some 
variable x is given by 

Since p — X^^Li M^s^, we find that 

^ = {<f>a\MM (5.23) 



and so, taking the derivative of eq.( 5.17 ), eq.( 5.2C| ) becomes 



a— 1 



Hence 



{ASf ^f:(Y^{c^a\MM^^) A.. (5.25) 

= l \a = l / 



For the experimental example given above, S = 0.106 ± 0.049. 



2. Linear Entropy 



The "linear entropy" is used to quantify the degree of mixture of a quantum state in an analytically convenient 
form, although unlike the von Neumann entropy it has no direct information theoretic implications. In a normalized 
form (defined so that its value lies between zero and one), the linear entropy for a two qubit system is defined by: 
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^ = i(l-Tr{p^}) 



(5.26) 



To calculate the error in this quantity, we need the foUowing partial derivative: 

r, 4 



dp 



O dpa 
a—1 



= -X^Pa(0a|M,|(/.,) 

--^r{pM.} 

„ 16 



(5.27) 



Hence the error in the Unear entropy is 



16 /g 16 \^ 



(5.28) 



For the example given in Sections 111 and IV, V = 0.037 ± 0.026. 



3. Concurrence, Entanglement of Formation and Tangle 

The concurrence, entanglement of formation and tangle are measures of the quantum-coherence properties of a 
mixed quantum state [ pof . For two qubits |^ , concurrence is defined as follows: consider the non-Hermitian matrix 
R = pYifp-Yj where the superscript T denotes transpose and the "spin flip matrix" E is defined by: 



/ -1 \ 
10 
10 
\ -1 j 



(5.29) 



Note that the definition of S depends on the basis chosen; we have assumed here the "computational basis" 
{\HH), \ HV), \ VH), \ VV)}. In what follows, it will be convenient to write R in the following form: 



- 1 

fl^U — l 



(5.30) 



where q^^i, — M^SMjS + MjySMjS. The left and right eigenstates and eigenvalues of the matrix R we shall denote 
by (Ca|, ICa) and r^, respectively, i.e.: 



{£,a\R^ra{£,a\ 
R\Ca)^ra\Ca). 



(5.31) 



t^The analysis in this subsection applies to the two qubit case only. Measures of entanglement for mixed n-qubit systems are 
a subject of on-going research: see, for example, [B0| for a recent survey. 
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We shall assume that these eigenstates are normalized in the usual manner for bi-orthogonal expansions, i.e. (^alCfe) = 
5a.b- Further wc shall assume that the eigenvalues are numbered in decreasing order, so that ri > r2 > rs > r4. The 
concurrence is then defined by the formula 



C = Max {0, Vn" - - - V^} 

= Max|o,^sgn('|-a) V^j, (5.32) 



k a=l ^ ^ ) 

where sgn(x) = 1 if a; > and sgn(a;) = — 1 if a; < 0. The tangle is given by T = and the Entanglement of 
Formation by 

E = ,{l±^y (5.33) 

where ft, (.t) = — .Tlog2.T— (1 — .T)log2(l — .t). Because is a monotonically increasing function, these three quantities 
are to some extent equivalent measures of the entanglement of a mixed state. 

To calculate the errors in these rather complicated functions, we must employ the perturbation theory for non- 
Hermitian matrices (see Appendix C for more details). We need to evaluate the following partial derivative, 

dC /3 \ 1 dra 



^sgn 



o=l 

4 16 



= X^X^sgn (I -aj ^^{^a\q^.,uS|,\Ca), (5.34) 



a=l /i=l 



where the function sgn(a;) is the sign of the quantity x: it takes the value 1 if a; > and —1 if a; < 0. Thus 
sgn (3/2 — a) is equal to +1 if a = 1 and —1 if a = 2, 3 or 4. Hence the error in the concurrence is 

16 / 4 16 /Q \ 1 \^ 

For our example the concurrence is 0.963 ± 0.018. 

Once we know the error in the concurrence, the errors in the tangle and the entanglement of formation can be 
found straightforwardly: 

AT = 2CAC (5.36) 



where h'{x) is the derivative of h(x). For our example the the tangle is 0.928±0.034 and the entanglement of formation 
is 0.947 ±0.025. 



VI. CONCLUSIONS 



In conclusion, wc have presented a technique for reconstructing density matrices of qubit systems, including a full 
error analysis. We have extended the latter through to calculation of quantities of interest in quantum information, 
such as the entropy and concurrence. Without loss of generality, we have used the example of polarization qubits of 
entangled photons, but we stress that these techniques can be adapted to any physical realization of qubits. 
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APPENDIX A: THE T-MATRICES 



One possible set of F-matrices are generators of SU{2) (g) SU{2), normalized so that the conditions given in eq,3.8 
are fulfilled. These matrices are: 



r. = 



ri3 = 



/o 1 o\ 

10 
1 

V 1 / 

/ 1 \ 

10 

10 

\ 1 / 

/O -i\ 



i 

\i / 

/O 1 \ 

10 

-1 

\0 -1 / 
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Ffi = 



ri4 = 
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1/ 



(A.l) 



As noted in the text, this is only one possible choice for these matrices, and the final results are independent of the 
choice. 



APPENDIX B: THE A/-MATRICES AND SOME OF THEIR PROPERTIES 

The M matrices, defined by eq.( |3.15 ), are as follows: 
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Ml = 



/ 2 -{l-i) -(1 + 1\ 

-{l + i) i 

-(l-i) -i 

\ 1 oy 

/o 1 \ 

i -{l + i) 

Q -i -(l-i) 
yi -(l-i) -(1 + 2 







2i 



Mr-, 



Mv = 



(l-i) 

-2i (l + i) 
\-(l-i) 



/ 

-(l-i) 

-(l + i) 
\-(l-i) 2 



-(1 + \ 



/ 



-(1 + \ 

2 


/ 



Mo 



M- 
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/ i \ 
-i 
i 

\ -i y 

/ i \ 
i 
-i 

V -i y 

/ 

-(l + i) 

-(l-i) 
\-(l-i) 2 



M 
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-(i + 0\ 


2 

/ 



Mi5 = 



/ -2i -(l-i)\ 

2i (l-i) 

(l+i) 

\-(i + «) y 



M2 = h 



M4 = i 



/ -(1-0 1 \ 

-(l + i) 2 i -(l + i) 

-i 

V 1 -(1 + y 

/ -(l + i) 1 \ 

i 

-(1 - -i 2 -(1 - i) 

\ 1 -(l + i) y 
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\-(l + i) 2i 



-(1 + \ 
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-(1-0 \ 



-2i 

y 



/ 1 \ 

0-10 
0-100 

\ 1 oy 



The form of these matrices is independent of the chosen set of matrices 



(B.l) 



used to convert the density matrix 



into a column vector. However the AIi, matrices do depend on the set of tomographic sta tes |'0 i/). 

There are some useful properties of these matrices which we will now derive. From eg. ( 3. 15 ), we have 



(B.2) 



From eq.( 3.12| ) we have (V'/ilf aIO/j) = ^/j.a, thus we obtain the result 



(B.3) 



If we denote the basis set for the four-dimensional Hilbert space by {|i) (i = 1, 2, 3, 4)}, then eq.(3.14) can be written 
as follows: 



(B.4) 



k,l V 
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Since cq.(B.4) is valid for arbitrary states p, we obtain the following relationship: 

Y,(fl\MA]){'^Ak)m.) = 5,k5fl. (B.5) 



Contracting eq.(B.5) over the indices («, j) we obtain: 

^Tr{M,}|V,)(7A,| =/, (B.6) 

where / is the identity operator for our four dimensional Hilbcrt space. 
A second relationship can be obtained by contracting eq.(|B.5|), viz: 



5](*|Af,|j)=<5„-, (B.7) 

or, in operator notation, 

^ iU = /. (B.8) 

APPENDIX C: PERTURBATION THEORY FOR NON-HERMITIAN MATRICES 

Whereas perturbation theory for Hermitian matrices is covered in most quantum mechanics text-books, the case of 
non-Hermitian matrices is less familiar, and so we will present it here. The problem is, given the eigenspectrum of a 
matrix _Ro [^ll' 

(.ia\R^ = r,i^U (C.l) 

i?0|Ca) = r.ICa), (C.2) 

where 

(?a|a>='5a,6 (C.3) 

we wish to find expressions for the eigenvalues r'^ and eigenstates (^^| and |Ca) of the perturbed matrix R' — Rq + SR. 

We start with the standard assumption of perturbation theory, i.e. that the perturbed quantities r^, (^^| and |Co) 
can be expressed as power series of some parameter A: 

r'a - -i"' + Ari^) + A^ri^) + . . . (C.4) 

IO-ie^>+A|e))+A2|C(2)> + ... (C.5) 

(CI-(e^l + A(e)| + A2(e^| + ... (C.6) 

Writing R' = Rq + XSR, and comparing terms of equal powers of A in the eigen equations, one obtains the following 
formulae: 

i?oie^)=ri")|Cr)) (C.7) 

(e^|i?o=ri")(e)| (C.8) 

(i?o-r(°)/)|e))=-(5i?-rW)|C(")) (C.9) 

(e^l(i?o-ri°)/) =-(e)|(<5i?-rW). (C.IO) 

Equations ( |C.7j) and ( |C.§| ) imply that, as might be expected, 

ld"^> = ICa) (C.ll) 

= (C.12) 

r^^^^r, (C.13) 
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Taking the inner product of eq.( |C.9D with and using the bi-orthogonal property eq.( p.3D , we obtain 

r^a^ = {USR\a- (C.14) 

This imphes that 



Thus, dividing both sides by some differential increment 5x and taking the hmit 8x — *■ 0, we obtain 

^-(f \ — \C) 



(C.15) 



(C.16) 



Using the completeness property of the eigenstates, |Cfc)(Cfc| — I a-nd the identity i?o — X^b ^^iCb) (Cf)| obtain 



the following formula 



b 



(C.17) 



Applying this to eq.(C.9) we obtain 



'{^b\SR\Ca) 



E 

b 

b^a 



n - Ta 



Kb). 



(C.18) 



Similarly, eqs.(C.lO) and (C.17) imply 



^{U6R\Cb) 



E 

6 

b^a 



(6 



(C.19) 
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FIG. 1. Tree diagram representing rmnibcr and type of measurements necessary for tomography. For a single qubit, 
the measurements {/to, Aii Aa, As} suffice to reconstruct the state, e.g., measurements of the horizontal, vertical, diagonal 
and right-circular polarization components, (H,V,D,R). For two qubits, 16 double-coincidence measurements are necessary 
({AoAo, AoAi • • • AsAsDj increasing to 64 three-coincidence measurements for three qubits ({AoAoAo, AoAoAi • • • AsAsAs})! and 
so on, as shown. 



"Black Box" 

source of photon pairs in 
arbitrary quantum states 



photon detector 
polarizer\ 




output 



polarizer- 
photon detector 

FIG. 2. Schematic illustration of the experimental arrangement. "QWP" stands for quarter waveplate, "HWP" for half 
waveplate; the angles of both pairs of wave plates can be set independently giving the experimenter four degrees of freedom 
with which to set the projection state. In the experiment, the polarizers were realized using polarizing prisms, arranged to 
transmit vertically polarized light. 
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Table 1 



TABLE 1: The tomographic analysis states used in our experiments. The number of coincidence counts measured 
in projections measurements provide a set of 16 data that allow the density matrix of the state of the two modes to 
be estimated. We have used the notation |D) = (|H) + |V)) /v^, |L) = (|H) + i|V)) /^/2 and |R) = (|H) - /v^. 
Note that, when the measurement are taken in the order given by the table, only one waveplate angle had to be 
changed between each measurement. 
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